%  SIGNALDISTORT Add distortion effects to a signal
%
% data = signaldistort(sig,
%                 ['BER',BER,]
%                 ['NOISELEVEL',NOISELEVEL,]
%                 ['NOISEWINDOW','windowtype']
%                 ['BURSTSTART',BURST START,'BURSTEND',BURST END],
%                 ['PLOTON'])
%
%   Note: Arguments in brackets [] are optional.
%
% SIG           -       Input signal
%
% BER           -       Bit Error Rate
%                       BER is a number in the range 0<=BER<=1.  It is the
%                       probability that a given bit will be flipped.
%
% NOISELEVEL    -       Amplitude of the background white noise
%
% NOISEWINDOW   -       White noise shaping window.  Options are:
%                           - 'none' : No shaping
%                           - 'rampUp' : noise ramped from 0% at start to
%                                        100% at end
%                           - 'rampDown': noise ramped from 100% at start
%                                         to 0% at end
%                           - 'notch' : noise ramped from 100% at start to
%                                       0% at middle to 100% at end
%                           - 'peak' : noise ramped from 0% at start to
%                                      100% at middle to 0% at end
%
% BURSTSTART    -       Start of noise 'burst' (high amplitude noise)
%                       region.  Burst start is a number in the range 
%                       0 <= N < 1, and is the percentage of the signal
%                       before the burst begins.  Burst_start must be less
%                       than burst_end.
%
% BURSTEND      -       End of noise 'burst' (high amplitude noise) region.
%                       Burst end is a number in the range 0 < N <= 1, and
%                       is the percentage of the signal before the burst
%                       begins.  Burst_end must be greater than
%                       burst_start.
%
% PLOTON        -       Plotting switch.  If 'ploton' is included,
%                       signalDistort generates a plot showing the signal
%                       before and after distortion.
%

function data = signalDistort(sig,varargin)
 
    %Init values
    BER = 0;
    noise_level = 0;
    window = 'NONE';
    burst_start = [];
    burst_end = [];
    ploton = false;

    %Grab the input args out of varargin  
    for i = 1:2:nargin-1
        switch upper(varargin{i})
            case 'BER'
                BER = varargin{i+1};
            case 'NOISELEVEL'
                noise_level = varargin{i+1};
            case 'NOISEWINDOW'
                window = upper(varargin{i+1});
            case 'BURSTSTART'
                burst_start = varargin{i+1};
            case 'BURSTEND'
                burst_end = varargin{i+1};
            case 'PLOTON'
                ploton = true; 
                i = i - 1; 
        end %end switch
    end %end for
    
    %Validate inputs and throw exceptions for bad inputs    
    %Check BER input
    if ~isnumeric(BER) || ischar(BER) %BER input isn't a number
        throw(MException('SignalDistorter:invalidInput',...
            sprintf('Error: Invalid input\nInput ''BER'' must have a numeric value')));
                
    elseif BER < 0 || BER > 1 %BER input is outside the valid range
        throw(MException('SignalDistorter:invalidInput',...
            sprintf('Error: Invalid input\nInput ''BER'' must lie in the range [0,1]')));
        
    %Check NOISELEVEL input
    elseif ~isnumeric(noise_level) || ischar(noise_level) %NOISELEVEL input value isn't a number
        throw(MException('SignalDistorter:invalidInput',...
            sprintf('Error: Invalid input\nInput ''noiselevel'' must have a numeric value')));
        
    %Check WINDOW
    elseif ~ismember(window,{'NONE','RAMPUP','RAMPDOWN','NOTCH','PEAK'})
        throw(MException('SignalDistorter:invalidInput',...
            sprintf('Error: Invalid input\nUnrecognized value for ''window''.  Run "signal_distort()" for more information')));
        
    %Check BURSTSTART and BURSTEND
    elseif (isempty(burst_start) && ~isempty(burst_end)) || (~isempty(burst_start) && isempty(burst_end))
        throw(MException('SignalDistorter:invalidInput',...
            sprintf('Error: Invalid input\nCannot define only 1 of ''BURSTSTART'' and ''BURSTEND''.  Instead, define both or neither')));
    
    elseif ~isempty(burst_start) && ~isempty(burst_end)
        
        if ~isnumeric(burst_start) || ~isnumeric(burst_end) ||...
            ischar(burst_start) || ischar(burst_end)
            throw(MException('SignalDistorter:invalidInput',...
                sprintf('Error: Invalid input\nInputs ''burst_start'' and ''burst_end'' must have numeric values')));

        elseif burst_start >= burst_end && burst_start >= 0 && burst_end > 0
            throw(MException('SignalDistorter:invalidInput',...
                sprintf('Error: Invalid input\nBurst start must be less than burst end')));

        elseif burst_start < 0 || burst_start >= 1
            throw(MException('SignalDistorter:invalidInput',...
                sprintf('Error: Invalid input\nInput ''burst_start'' must lie in the range (0,1]')));

        elseif burst_end <= 0 || burst_end > 1
            throw(MException('SignalDistorter:invalidInput',...
                sprintf('Error: Invalid input\nInput ''burst_end'' must lie in the range [0,1)')));
        end %end if
    end %end if
    
    %Load signal and convert to lin-16 if necessary
    data = sig;
    
    %%%IMPORTANT - Store original data in "origData".  This allows the
    %%%various impairments to be added at the end instead of being applied
    %%%throughout.
    origData = data;
    
    %If plotting is enabled, save the "before" signal and spectrogram
    if ploton
        beforeData = data;        
    end %end if

    %Apply bit errors
    if BER > 0
        nbits = 8;
        
        %Create random array of bit flip probabilities
        R1 = rand(size(data));
        
        %Threshold based on BER        
        R1(R1 <= BER) = 0;
        R1(R1 > BER) = 1;
        R1 = ~R1;
        
        %Create another random array of bit flip locations and signs
        R2 = floor(rand(size(data)) * nbits) .* sign(randn(size(data)));
        
        %Combine the two arrays into a bitmask
        bitmask = R1 .* 2.^R2;
        
        data = data + bitmask;
    end %end if 
    
    %Apply bursted noise
    if ~isempty(burst_start) && ~isempty(burst_end)
        windim = [floor(burst_start * length(data))+1 floor(burst_end*length(data))+1];
        bwin = zeros(size(data));
        bwin(windim(1):windim(2)) = randn(windim(2)-windim(1)+1,1) .* ...
            hamming(windim(2)-windim(1)+1) * 0.2;
        data = data + bwin;        
    end %end if
    
    %Create background noise, apply window to noise level, and apply noise
    if noise_level > 0
        noise = randn(size(data));
        noise = noise_level * (noise / max(noise));
        switch window
            case 'RAMPUP'
                noise = noise .* linspace(0,1,length(noise))';
            case 'RAMPDOWN'
                noise = noise .* linspace(1,0,length(noise))';
            case 'NOTCH'
                noise = noise .* [linspace(1,0,floor(length(noise)/2)) ...
                                  linspace(0,1,ceil(length(noise)/2))]';
            case 'PEAK'
                noise = noise .* [linspace(0,1,floor(length(noise)/2)) ...
                                  linspace(1,0,ceil(length(noise)/2))]';
        end %end switch
        
        data = data + noise;
    end %end if
    
    %Renormalize to [-1,1]
    data = data / max(abs(data));

     %If plotting, generate plots
    if ploton
        figure('Name','Signal And Distortion Effects','Numbertitle','off');
        
        %Plot time domain of "before"
        subplot(211)
        plot(beforeData)
        xl = xlim;
        ylabel('Amplitude');
        grid on
        
        %Plot time domain of "after"
        subplot(212)
        plot(data)
        xl = xlim;
        ylabel('Amplitude');
        grid on
        
        %"Maximize" the figure (not true Windows maximized, just full
        % screen)
        units=get(gcf,'units');
        set(gcf,'units','normalized','outerposition',[0 0 1 1]);
        set(gcf,'units',units);

    end %end if
    
end %end function
